Modeling the Volatility of U.S Bond Yields

The prompt and base code for this project can be found on Datacamp. The basis of this is to explore the volatility of Zero Coupon Bond Yields (from year maturities 1-30) and see how that has changed since the 1960’s. Below, I loaded the necessary librares and did some data exploration of the time series data obtained from Quandl. This was obtained making an API call.

Necessary packages

defaultW <- getOption("warn") 

options(warn = -1) 

library(dygraphs)
library(Quandl)
library(dplyr)
library(rugarch)
library(plotly)
library(RColorBrewer)

options(warn = defaultW)

#get dataset, filter for coronavirus dates
df <- Quandl("FED/SVENY", api_key="15cxBvvbCzucYDswsDfJ")
covid_dates <- subset(df, Date > '2020-01-01' & Date < '2020-08-20')

str(df)

'data.frame':   14755 obs. of  31 variables:
 $ Date   : Date, format: "2020-08-14" "2020-08-13" "2020-08-12" "2020-08-11" ...
 $ SVENY01: num  0.157 0.162 0.162 0.159 0.146 ...
 $ SVENY02: num  0.157 0.171 0.163 0.155 0.135 ...
 $ SVENY03: num  0.197 0.215 0.202 0.187 0.161 ...
 $ SVENY04: num  0.26 0.278 0.261 0.24 0.209 ...
 $ SVENY05: num  0.334 0.351 0.331 0.304 0.268 ...
 $ SVENY06: num  0.412 0.426 0.404 0.373 0.332 ...
 $ SVENY07: num  0.491 0.502 0.478 0.443 0.398 ...
 $ SVENY08: num  0.568 0.577 0.551 0.512 0.463 ...
 $ SVENY09: num  0.643 0.649 0.621 0.58 0.527 ...
 $ SVENY10: num  0.715 0.717 0.688 0.644 0.589 ...
 $ SVENY11: num  0.784 0.783 0.752 0.706 0.649 ...
 $ SVENY12: num  0.849 0.846 0.813 0.766 0.706 ...
 $ SVENY13: num  0.911 0.906 0.872 0.823 0.762 ...
 $ SVENY14: num  0.971 0.964 0.927 0.878 0.815 ...
 $ SVENY15: num  1.028 1.019 0.98 0.93 0.866 ...
 $ SVENY16: num  1.082 1.071 1.031 0.98 0.915 ...
 $ SVENY17: num  1.134 1.121 1.08 1.028 0.962 ...
 $ SVENY18: num  1.18 1.17 1.13 1.07 1.01 ...
 $ SVENY19: num  1.23 1.22 1.17 1.12 1.05 ...
 $ SVENY20: num  1.28 1.26 1.21 1.16 1.09 ...
 $ SVENY21: num  1.32 1.3 1.25 1.2 1.13 ...
 $ SVENY22: num  1.36 1.34 1.29 1.24 1.17 ...
 $ SVENY23: num  1.4 1.38 1.33 1.28 1.21 ...
 $ SVENY24: num  1.44 1.42 1.37 1.31 1.25 ...
 $ SVENY25: num  1.48 1.45 1.4 1.35 1.28 ...
 $ SVENY26: num  1.52 1.49 1.43 1.38 1.31 ...
 $ SVENY27: num  1.55 1.52 1.46 1.41 1.34 ...
 $ SVENY28: num  1.58 1.55 1.5 1.44 1.38 ...
 $ SVENY29: num  1.61 1.58 1.52 1.47 1.4 ...
 $ SVENY30: num  1.64 1.61 1.55 1.5 1.43 ...
 - attr(*, "freq")= chr "daily"

Plotting daily estimates for yall zero coupon yields for 2020

dfasxts <- as.xts(x = df[, -1], order.by = df$Date)
covid_dates <- as.xts(x = covid_dates[, -1], order.by = covid_dates$Date)
dygraph(covid_dates, main = "All Zero Coupon Yields (1-30) 2020", ylab = "Value") %>%
            dyAxis('x', axisLabelFontSize = 12) %>%
            dyRangeSelector()

df$Date <- as.Date(df$Date)
df$year <- format(as.Date(df$Date, format="%m/%d/%Y"),"%Y")
df <- select(df, -Date)
df <- na.omit(df)

df <- df %>%
  group_by(year) %>%
  summarise_all(mean)

  
SVENY01 <- df$SVENY01
SVENY10 <- df$SVENY10
SVENY30 <- df$SVENY30


data <- data.frame(df, SVENY01, SVENY10, SVENY30)
data <- select(data, year, SVENY01, SVENY10, SVENY30)
fig <- plot_ly(data, x = data$year, y = ~SVENY01, name = 'SVENY01', type = 'scatter', mode = 'lines + markers') 
fig <- fig %>% add_trace(y = ~SVENY10, name = 'SVENY10', mode = 'lines + markers') 
fig <- fig %>% add_trace(y = ~SVENY30, name = 'SVENY30', mode = 'lines + markers')
fig <- fig %>% layout(title = "",
         xaxis = list(title = "year"),
         yaxis = list (title = "value"))

1, 10, and 30 year maturity bond yield averages across time

fig

Bond yield heat map

row.names(df) <- df$year

Setting row names on a tibble is deprecated.

df <- select(df, -year)
df_matrix <- data.matrix(df)

df_heatmap <- heatmap(df_matrix, Rowv=NA, Colv=NA, col = brewer.pal(6, "Greys"), scale="column", margins=c(4,1), main = "Heatmap")

# plotting the evaluation of bond yields
library(viridisLite)
yields <- dfasxts
plot.type <- "single"
plot.palette <- magma(n = 30)
asset.names <- colnames(dfasxts)
plot.zoo(x = dfasxts, plot.type = "single", col = plot.palette, ylab = "", xlab = "")
legend(x = "topleft", legend = asset.names, col = plot.palette, cex = 0.45, lwd = 3)

dfasxts_d <- diff(dfasxts)

plot.zoo(x = dfasxts_d, plot.type = "multiple", ylim = c(-0.5, 0.5), cex.axis = 0.7, ylab = 1:30, col = plot.palette, main = "", xlab = "")

dfasxts <- dfasxts_d["2000/",]
x_1 <- dfasxts[,"SVENY01"]
x_20 <- dfasxts[, "SVENY20"]

# Plot the autocorrelations of the yield changes)
par(mar=c(5.1, 4.1, 4.1, 2.1))
par(mfrow=c(2,2))
acf_1 <- acf(x_1)
acf_20 <- acf(x_20)

# Plot the autocorrelations of the absolute changes of yields
acf_abs_1 <- acf(abs(x_1))
acf_abs_20 <- acf(abs(x_20))

NA
NA

spec <- ugarchspec(distribution.model = "sstd")


fit_1 <- ugarchfit(x_1, spec = spec)


vol_1 <- sigma(fit_1)
res_1 <- scale(residuals(fit_1, standardize = TRUE)) * sd(x_1) + mean(x_1)


merge_1 <- merge.xts(x_1, vol_1, res_1)
plot.zoo(merge_1, xlab = "Year")

####################

fit_20 <- ugarchfit(x_20, spec = spec)


vol_20 <- sigma(fit_20)
res_20 <- scale(residuals(fit_20, standardize = TRUE)) * sd(x_20) + mean(x_20)


merge_20 <- merge.xts(x_20, vol_20, res_20)
plot.zoo(merge_20, xlab = "Year")

par(mar=c(5.1, 4.1, 4.1, 2.1))
par(mfrow=c(2,1))
hist(res_1)
hist(res_20)

ugarchspec()

*---------------------------------*
*       GARCH Model Spec          *
*---------------------------------*

Conditional Variance Dynamics   
------------------------------------
GARCH Model     : sGARCH(1,1)
Variance Targeting  : FALSE 

Conditional Mean Dynamics
------------------------------------
Mean Model      : ARFIMA(1,0,1)
Include Mean        : TRUE 
GARCH-in-Mean       : FALSE 

Conditional Distribution
------------------------------------
Distribution    :  norm 
Includes Skew   :  FALSE 
Includes Shape  :  FALSE 
Includes Lambda :  FALSE 

density_x_1 <- density(x_1)
density_res_1 <- density(res_1)


plot(density_x_1)
lines(density_res_1, col = "red")

norm_dist <- dnorm(seq(-0.4, 0.4, by = .01), mean = mean(x_1), sd = sd(x_1))
lines(seq(-0.4, 0.4, by = .01), 
      norm_dist, 
      col = "darkblue"
     )

# Add legend
legend <- c("Before GARCH", "After GARCH", "Normal distribution")
legend("topleft", legend = legend, 
       col = c("black", "red", "darkblue"), lty=c(1,1))

distribution <- qnorm


qqnorm(x_1, ylim = c(-0.5, 0.5))
qqline(x_1, distribution = distribution, col = "darkgreen")

par(new=TRUE)
qqnorm(res_1 * 0.614256270265139, col = "red", ylim = c(-0.5, 0.5))

qqline(res_1 * 0.614256270265139, distribution = distribution, col = "darkgreen")
legend("topleft", c("Before GARCH", "After GARCH"), col = c("black", "red"), pch=c(1,1))

Conclusion

Garch modeling brought the residuals closer to normal distribution. Year 1 yield deviates more from a normally distributed white noise process than a 20 year yield.